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Abstract. An analytical scheme and a numerical method in 
order to study the effects of general relativity on the viscos- 
ity driven secular bar mode instability of rapidly rotating stars 
are presented. The approach consists in perturbing an axisym- 
metric and stationary configuration and studying its evolution 
by constructing a series of triaxial quasi-equilibrium configu- 
rations. These are obtained by solution of an approximate set 
of field equations where only the dominant non-axisymmetric 
terms are taken into account. The progress with respect to 
our former investigation consists in a higher relativistic order 
of the non-axisymmetric terms included into the computation, 
namely the fully three-dimensional treatment of the vector part 
of the space-time metric tensor as opposed to the scalar part, 
solely, in the former case. The scheme is applied to rotating 
stars built on a polytropic equation of state and compared to 
our previous results. The 3D-vector part turns out to inhibit the 
symmetry breaking efficiently. Nevertheless, the bar mode in- 
stability is still possible for an astrophysically relevant mass of 
M ns = 1.4 Mq when a stiff polytropic equation of state with an 
adiabatic index of 7 = 2.5 is employed. Triaxial neutron stars 
may be efficient emitters of gravitational waves and are thus po- 
tentially interesting sources for the forthcoming laser interfero- 
metric gravitational wave detectors such as LIGO, VIRGO and 
GEO600. From a numerical point of view, the solution of the 
three-dimensional minimal-distortion shift vector equation in 
spherical coordinates is an important achievement of our code. 

Key words: gravitation - instabilities - relativity - stars: inte- 
riors - stars: rotation 



1. Introduction 

Large scale laser interferometric gravitational detectors such 
as LIGO, VIRGO and GEO600 are supposed to become op- 
erative in a few years, and among the possible astrophysical 
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scenarios that of a rapidly rotating neutron star, deviating pro- 
gressively from the initially axisymmetric and stationary con- 
figuration towards a triaxial one, is of particular interest. In this 
case, a "bar shaped" neutron star could become an efficient 
source of continuous wave gravitational radiation. Rapidly ro- 
tating stars can break their symmetry, if the ratio T/|W| of 
kinetic and potential energy exceeds some critical value (see 
Schutz 1987 for a review) and the (/ = 2, m = 2) bar mode 
gets unstable. This can be the case for a newborn neutron star, 
having acquired a sufficiently large amount of rotational kinetic 
energy during its formation. The matter viscosity being low, it 
may be sensitive to the Chandrasekhar-Friedman-Schutz insta- 
bility which is related to gravitational radiation reaction (Wag- 
oner 1984; Lai & Shapiro 1995). An alternative mechanism 
is the viscosity driven secular instability which is eventually 
operational in old neutron stars in close binary systems, be- 
ing spun up by matter accretion from a nearby companion to 
some critical value of T/|W| (Ipser & Managan 1984; Bonaz- 
zola et al. 1996). As concerns rigidly rotating, homogeneous 
and self-gravitating fluid bodies, Newtonian theory tells that at 
moderate angular velocity, the fluid body is shaped like some 
axisymmetric Maclaurin spheroid. At T/|W| = 0.2738, the 
Maclaurin spheroids become dynamically unstable. However, 
at T/\W I = 0.1375 exists a bifurcation point towards two fam- 
ilies of triaxial configurations via the above introduced secular 
instabilities. In the case of the viscosity driven instability, ki- 
netic energy is dissipated whereas angular momentum is con- 
served. As a consequence, the fluid body evolves along a se- 
quence close to some Riemann S ellipsoids towards a Jacobi 
ellipsoid which is the configuration with the lowest rotational 
kinetic energy for the given angular momentum and mass of 
the initial Maclaurin spheroid (Press & Teukolsky 1973). The 
transition towards the Jacobi ellipsoid occurs on the associated 
viscous time-scale which is much longer than the dynamical 
one. For this reason, it is called a secular instability. In the 
final state, viscous dissipation has ceased, and the fluid body 
rotates rigidly again about its smallest axis in an inertial frame. 
It is this viscosity driven secular instability which we are con- 
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cerned with in the following. In the Newtonian compressible 
case, it was shown by Jeans (1919, 1928), that for a polytropic 
equation of state 7 2.2 is needed to reach the bifurcation 
point. If the equation of state is softer, the critical angular ve- 
locity O cr it is greater than the mass shedding limit £!k, and 
the bifurcation point is inaccessible. James (1964) has given a 
refined value for the critical adiabatic index of 7 crit = 2.238, 
a value which has been recently confirmed by Bonazzola et al. 
(1996) and, independently, by Skinner & Lindblom (1996) who 
have found a corresponding value of 7 cr j t = 2.237. The result of 
Bonazzola et al. (1996) was a by-product of the investigation of 
the spontaneous symmetry breaking of fully relativistic rapidly 
rotating stars, presented in the same paper. Until then, all stud- 
ies of the bifurcation point had been carried out for Newto- 
nian configurations (Jeans 1928; James 1964; Ipser & Man- 
agan 1981; Hachisu & Eriguchi 1982; Skinner & Lindblom 
1996) or at the 1-PN level of a post-Newtonian analysis, at most 
(Chandrasekhar 1967; Tsirulev & Tsvetkov 1982). In view of 
the highly relativistic character of realistic neutron star mod- 
els, this approach is rather motivated by technical simplifica- 
tion than by physical reasoning. Our approach consists in per- 
turbing the (7 = 2, to = 2) bar mode of an "exact" axisymmet- 
ric configuration and studying the growth or the decay of the 
applied perturbation. We take into account only the dominant 
non-axisymmetric terms up to a certain order. In Bonazzola et 
al. (1996), the level of approximation corresponded to an ex- 
pansion of the Einstein equations of order 0-PN of the three- 
dimensional terms. In the present paper, we have included the 
fully three-dimensional treatment of the shift vector N l , rais- 
ing the approximation level to order 1/2-PN. We recall the basic 
assumption of rigid rotation in our approximation (see Bonaz- 
zola et al. (1996) for a discussion of the astrophysical context), 
and the negligibility of gravitational radiation. Approximation 
schemes which exploit the near equilibrium of the gravitational 
fields, have gained increasing interest in the context of the bi- 
nary coalescence problem. Wilson & Mathews (1989) were the 
first ones to propose the solution of dynamical problems in gen- 
eral relativity by integration of a reduced set of field equations, 
which basically corresponds to the successive solution of initial 
value problems for the gravitational fields (York 1979) whereas 
the hydrodynamic part is calculated exactly under the action of 
the momentary gravitational fields. It has been applied in the 
following to study the late stages of an inspiralling neutron star 
binary system (Wilson & Mathews 1995; Wilson et al. 1996). 
Cook et al. (1996), being interested in the initial value problem 
of coalescing binaries, have investigated the validity of Wil- 
son's approach by application to stationary and axisymmetric 
neutron star models which, by construction, are perfect equi- 
librium systems. They further modified Wilson's approach in 
a way to exploit the near equilibrium of the matter fields as 
well, which allows to reduce the matter evolution equations to 
an algebraic first integral of motion. The neutron star models 
computed by means of the simplified field equations coincide 
with those based on the full set of Einstein equations within a 
few percent. One of the essential points of the Wilson approach 
is to adopt a conformally flat three-metric in order to simplify 



the field equations. This "conformally flat condition" is being 
justified by the negligible impact of the radiation part of the 
gravitational fields on the matter motion. Indeed, for problems 
in spherical symmetry where no radiative component exists, it 
is possible to adopt isotropic coordinates without any analytic 
approximation. However, it has to be pointed out, that even for 
stationary and axisymmetric systems where no radiative part is 
present either, isotropic coordinates are not sufficient to give a 
complete description of the problem. A further degree of free- 
dom of the gravitational field is required which is subject to the 
dynamical Einstein equations. The other gravitational fields are 
then entirely amenable to determination by means of elliptic 
equations derived from (1) the Hamiltonian constraint equation 
for the conformal factor, (2) the momentum constraint equa- 
tion for the shift vector, and (3) the additional maximal slicing 
condition which determines the lapse function. Our model im- 
proves on Wilson's scheme in so far as our approximate equa- 
tions reduce to the exact Einstein equations, as long as the con- 
figuration remains stationary and axisymmetric, so that only 
the three-dimensional perturbation is treated approximately. 

The paper is organized as follows: In Sect. 2 we present the 
basic assumptions of our model and derive the field and matter 
equations. Sect. 3 introduces the numerical code, and in Sect. 4 
we present the results of the improved scheme for a polytropic 
equation of state, including the fully three-dimensional solu- 
tion of the shift vector equation, as well as a comparison with 
the results of our previous study. In Sect. 5 we will draw some 
concluding remarks. 

2. Theoretical model 

2.1. Basic assumptions 

The general space-time line element in terms of the quantities 
of the (3+l)-formalism of general relativity (see e.g. York 1979 
for an introduction) is given by 

Ax^&x" = -N 2 dt 2 + h t] (dx l -ArMi)(dx J -A"dt) ,(1) 

with the lapse function N, the shift vector N l , and hij, the 
metric tensor induced in the spatial hypersurfaces E t . Before 
the symmetry breaking, the space-time associated with the ro- 
tating star is stationary and axisymmetric. We briefly recall the 
main conclusions for the case where the star matter is assumed 
to be constituted by a perfect fluid, the stress energy tensor hav- 
ing the form 

T =(e+p)u®u + pg . (2) 

The involved quantities are the fluid proper energy density e, 
the fluid pressure p, the fluid four-velocity u, and the space- 
time metric tensor g. Two Killing vector fields k and m are 
linked to the space-time symmetries where k is time-like at 
least far from the star and I space-like, its orbits being closed 
curves. In the case of rigid rotation, the space-time is circular 
and the two-surfaces orthogonal to both k and m are globally 
integrable (Carter 1973). The coordinates t and <p are associated 
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with the both Killing vector fields k~dt and m = d ( j > whereas 
the remaining coordinates r and 8 can be chosen arbitrarily. 
The standard coordinates for stationary axisymmetric systems 
are quasi-isotropic coordinates where a conformally flat met- 
ric in the (r, 8) -coordinate planes is adopted. The general line 
element (1) specified to these coordinates reads 



g„ v dx»dx v = -N 2 dt 2 + B 2 r 2 sin 2 8 (d<f>- N^dt) 2 



(3) 



+A 2 (dr 2 



l d8 2 



The shift vector has only one non-vanishing component N^, 
which represents the dragging of inertial frames by the rotating 
star. For this coordinate choice, Bonazzola et al. (1993) have 
exhibited a compact set of elliptic equations for the metric po- 
tentials N, and A, B, the latter being defined by A = AN, 
B = BN . Let us also remind that quasi-isotropic coordinates 
satisfy both the minimal distortion coordinate condition, intro- 
duced by Smarr & York (1978), and the maximal slicing condi- 
tion K = 0, where K is the scalar of extrinsic curvature. When 
the star deviates from axisymmetry, it is no more stationary ei- 
ther, as gravitational radiation carries away energy and angular 
momentum. However, at the very beginning of the symmetry 
breaking, this deviation is very small, and, consequently, the 
losses due to gravitational radiation are negligible. Under this 
assumption, and for rigid rotation, there exists a Killing vec- 
tor field I which is proportional to the fluid velocity u (Carter 
1979), 



u = XI 



(4) 



where A is a strictly positive scalar function. In the stationary 
and axisymmetric case, the Killing vector I is given by 



I = k + f!m . 



(5) 



The constant ft is the angular velocity defined as Q = vr /vf. 
In the non-axisymmetric case, we assume (1) the existence of 
a vector field k which is time-like at least far from the star, (2) 
a vector field m which is space-like everywhere, (3) a constant 
Q, such that I defined by (5) is a Killing vector field and (4) the 
fluid velocity u is proportional to I. The Killing vector field I 
is associated with the persisting helical symmetry of the space- 
time generated by the non-axisymmetric body which appears 
still static in the corotating frame. 

2.2. Matter equations 

Based on the assumptions made in the previous section, namely 
that (1) the star matter is composed of a single constituent per- 
fect fluid, and (2) the star rotates rigidly, it is possible to derive 
a simple first integral of motion following the procedure out- 
lined in Bonazzola et al. (1996). We first introduce the family 
of Eulerian observers Co, whose four-velocity coincides with 
the future directed unit vector field n orthogonal to the space- 
like hypersurfaces S t . The Lorentz factor r between these lo- 
cal rest observers C and the fluid comoving observers 0\ is 
given by 



r 



-n ■ u . 



(6) 



With the baryon chemical potential p and the mean baryon 
mass tub, the log-enthalpy H is defined as 



\mflC 2 



(7) 



which is the relativistic generalization of the Newtonian spe- 
cific enthalpy h. Introducing v = In N, we recover the first in- 
tegral of motion 



H + v — In r = const , 



(8) 



already familiar from the axisymmetric and stationary case. 
Note, however, that in the present case all quantities are func- 
tions of (r, 8, ip) where ip = <fi — VL t is the azimuthal angular 
variable in the corotating frame . At the Newtonian limit, we 
have H -> h, v -> U, -lnT -» -l/2Vt 2 p 2 , and (8) ap- 
proaches the classical first integral of motion where U is the 
Newtonian potential and p the distance from the rotation axis. 

2.3. Field equations 

As announced in Sect. 1, we apply an approximate set of field 
equations derived under the assumptions that (1) the helical 
symmetry of space-time is conserved after deviation from the 
axisymmetric and stationary configuration, and (2) gravita- 
tional radiation is negligible. 

In addition, we only retain the dominant non-axisymmetric 
contributions in the field equations up to order 1/2-PN, their 
leading relativistic order being less or equal than a 3 / 2 where 



(9) 



is the post-Newtonian expansion parameter. At this level of ap- 
proximation, the lapse function N(r, 8, ip) and the shift vector 
N % (r,8,iji) have to be considered as three-dimensional quan- 
tities. The components N r (r,9,tp) and N e (r,8,ip) are gen- 
uinely three-dimensional contributions and are absent at the 
previous approximation level 0-PN. Corrections of higher rel- 
ativistic order to the metric tensor are included for the diag- 
onal components via the axisymmetric potentials A(r, 8) and 
B(r 7 8), their sources being essentially dominated by the fluid 
pressure whereas the extra-diagonal terms are again generically 
non-axisymmetric quantities and hence are neglected. Accord- 
ingly, the spatial metric tensor reads 



'A 2 



(hij) = 



N 2 



(10) 







B 2 r 2 sin 2 i 



The choice of A^ 2 as conformal factor of the (nearly flat) con- 
formal three-metric has proven to be particularly advantageous, 
as has been exposed by Bonazzola et al. (1993). It isolates 
the lapse function as the predominant part of the gravitational 
fields, which is underlined by considering the weak field limit 
of the corresponding space-time line element, given by 



g^dx^dx" = -(l-2v)dt 2 + (l + 2v) 

x(dr 2 + r 2 d8 2 + r 2 sin 2 8d(j> 2 ) , 



(11) 
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which conducts to the Newtonian equation for the gravitational 
potential v. 

Having specified the space-time line element for the per- 
turbed neutron star models, we can proceed to derive the gov- 
erning field equations. The temporal evolution of hij is deter- 
mined by 



dthij + N i{j + N j{l + 2NK tJ = , 



(12) 



with the tensor of extrinsic curvature Kij . For h = | hi j \ \ f% j | ~ 1 , 
the determinant of the spatial metric tensor, normalized by that 
of flat space spherical coordinates = r 4 sin 2 9, it follows 
immediately the relation 



d t h + 2(N l ll + NK)h = 



(13) 



where K = K j j denotes the trace of K iy Equations (12) and 
(13) enable us to derive the evolution equation of the conformal 
metric tensor ha = h^ 1 ^ ha 



h^dtlh^ 3 h^] + [Ni U + N Mi - f(A^) h^] (14) 
+2N [Kij-\Khij] = . 



The lapse function N, the only three dimensional quantity in- 
volved in the product /i -1 / 3 hij, cancels out in this term. There- 
fore, the temporal derivative of hij vanishes identically. If we 
further impose the maximal slicing condition K = 0, we can 
determine Kij from the metric potentials according to 



K ll = -[N l y + N 3ll -l(N l ll )h ll ]/2N. 



(15) 



Furthermore, the maximal slicing condition yields an elliptic 
equation for the lapse function N 



- N[Aw(E+S) + K kl K lk ] = . 



(16) 



Here E stands for the total energy density and S for the trace 
of the stress tensor S l j, all of them measured by the Eulerian 
observer O - 

Inserting (15) into the momentum constraint equation 



K ij \j -8irJ l = 



(17) 



leads immediately to the maximal slicing-minimal distortion 
shift vector equation 

Ni]j \j + \{N J \ j )\ i + R l j N 3 + 2K^N {j + WirNJ 1 = , (18) 

introduced by Smarr & York (1978) where J 1 denotes the 
momentum density vector. Indeed, the York minimal distor- 
tion gauge condition [^(/i 1 / 3 h l i)]\j = is trivially fulfilled, 
since already the interior of the square brackets equals 0. Note 
that any coordinate system, whose conformal metric tensor is 
time-independent, automatically satisfies the minimal distor- 
tion gauge condition. This is notably the case for isotropic 
coordinates, but also for our choice of quasi-isotropic coordi- 
nates where the conformal metric is not that of flat space like in 
the Wilson scheme, but time-independent as well. As a conse- 
quence, our approximation of keeping the original form of the 



spatial metric except for the lapse function N, being treated as 
a three-dimensional quantity now, ensures the coordinates to 
remain maximal slicing-minimal distortion coordinates in the 
three dimensional case after deviation from the initial station- 
ary and axisymmetric configuration. 

The explicit field equations for our particular choice ( 1 0) of 
the spatial metric tensor are then derived after introduction of 
the auxiliary variables 



ln^, /3 = lnB, and G = B r sin 9 



(19) 



following the procedure outlined in Bonazzola et al. (1993). 
We obtain the following elliptic equation for the logarithm v of 
the lapse function N 



A 3 ^ = tor^E+S) 
B 2 



(20) 



+ — r 2 sm 2 9 {dN 4 ') 2 -dud/3, 



where A 3 denotes the three-dimensional flat space scalar 
Laplacian with respect to the coordinates (r, 9, ip) of the coro- 
tating frame 



A, = 



d 2 



2d Id 2 



1 



d 



dr 2 r dr r 2 d0 2 r 2 tan 9 d9 
1 d 2 

+ - 



(21) 



r 2 sin 2 9 dip 2 
and where the following abridged notation 

ft,^=— ^ + - — ^ 
dr dr r 2 d9 d9 



(22) 



is used. We define the pseudo-physical components of the shift 
vector N l via the following relations 



(23) 



The shift vector equation (18) associated with this particular 
frame yields 



(24) 



with specified by (30) below. We further introduce the ex- 
plicit expressions of the involved derivative operators, associ- 
ated with the standard orthonormal frame of flat space spheri- 
cal coordinates. The pseudo-physical components of the three- 
dimensional flat space vector Laplacian AV are specified as 



(Avy = A 3 V f - 



2V r 

— - 

dV^ 



d_ _L_ 

d9 + tan 9 



V" 



(25) 



r 2 sin 9 dip 



(Ay)' 



r 2 sin 2 ( 



+ ■ 



2 dV f 



d9 



(26) 



2 



dV* 



r 2 sin 9 tan 9 dip 
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(Avy 



dV f 



r 2 sin 2 9 r 2 sin 9 dtp 
2 dV S 



-T7T- (27) 



r 2 sin 9 tan 9 dtp 

We further need to compute the covariant divergence VqV^ 
which reads 

„ TrS dV f 2V f 1 dV s V § 1 8V+ mm 

or r r d9 rt&nO rsmO dtp 

Finally, the gradient of a scalar potential U is computed accord- 
ing to 

or r 09 rsmv dip 

The pseudo-physical components of the actual source, com- 
puted by means of (18), read 



A6nNJ r , S tt 



= -16tt^ 



NA 2 J 



B 2 r sin 9 



-16nN— , and 



-rsin0<9A^<9(3/?-4i/) . 



(30) 



Equations (20) and (24) together with (30) constitute the 3D- 
part of our field equations. The remaining gravitational poten- 
tials A and B are computed by means of the dynamical Einstein 
and the Hamiltonian constraint equations after integration over 
ip, which conducts to the original equations derived in Bonaz- 
zola et al. (1993). They are genuine 2D-equations, intimately 
related to the axisymmetry and stationarity of the initial con- 
figuration (see Gourgoulhon & Bonazzola (1993) for a geo- 
metrically motivated derivation of the (3+l)-equations for this 
case). For the potentials a and G we then have 

A2& = + ^ r 2 sin 2 {dN*f (0„)A ,(31) 



A 2 G= Uir^-r sin6 (S r r +S e g) 



(32) 



where A 2 stands for the two-dimensional flat space scalar 
Laplacian 



A 2 EE 



d 2 



+ 



Id Id 2 



dr 2 r dr 



+ 



' d9 2 



(33) 



and where (A) denotes the average of A with respect to the 
angular variable ip. This ensures the consistency of the actual 
sources of (31) and (32) with the axisymmetry of a and G. 

To complete the analytic description, we add the matter re- 
lated quantities for a perfect fluid whose stress-energy tensor 
has been defined by (2), expressed in terms of variables of the 
(3+l)-formalism. With the Lorentz factor J 1 , defined by (6), 
the "physical" fluid velocity Ua with respect to the Eulerian 
observer O along the a-th coordinate line is given by 



(34) 



where ea is the corresponding spatial unit vector. For our coor- 
dinates, the components then read 



-_rAT*,and 



(35) 



and the normalization condition u ■ u = — 1 on the fluid four- 
velocity u yields 



r = (i - u 2 )- 1 ' 2 . 



(36) 



The matter related variables f, E, J, and S z j, specified to our 
coordinate system, take the approximate form 

r = (1 - U 2 )- 1 ' 2 , (37) 
E = r 2 (e+p)-p, (38) 

J lP = (E+p)^r 2 sm 2 9(n-N*) , (39) 

S r r =p, S\=p, S% = (E+p)U 2 +p, (40) 

whereas the remaining components equal 0. By combining (30) 
and (39) one obtains the final form of S?, 

S f = , S § = , and (41) 
A 2 

S 10 = 16tt(E+ P )— rsm9(n-N' t >) 
— r sin 

where the latter is exactly the same expression as the one pre- 
sented in Bonazzola et al. (1993) except that the lapse function 
TV, the shift vector component and the matter term (E+p) 
are allowed to depend on ip now. 

Let us finally mention that our analytic scheme yields the 
exact solution to the general field and matter equations for two 
limiting cases: (1) at the Newtonian limit for an arbitrary devi- 
ation from axisymmetry, and (2) in the axisymmetric case up 
to arbitrary relativistic order. 

2.4. Stability of an axisymmetric configuration 

At this point, we can summarize our analytical approach. The 
elliptic field equations (20), (24), (31), and (32), completed by 
the first integral equation (8), fully determine an axisymmet- 
ric and stationary equilibrium model, having specified e.g. the 
central value of the log-enthalpy H c and the angular velocity O 
for some particular equation of state. The above problem can 
be formulated as a fixed point problem in some appropriate 
functional Banach space. Under reasonable physical assump- 
tions, the induced mapping C is contractive, thus a unique so- 
lution exists and the deviation of the sequence members from 
the fixed point is bounded by some decaying exponential func- 
tion. We refer to Schaudt & Pfister (1996) for a recent proof of 
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this statement, though, at present, restricted to weakly relativis- 
tic configurations such as white dwarfs. The solution scheme 
consists in solving the three-dimensional field and matter equa- 
tions iteratively where as initial guess a spherical, static matter 
distribution with a parabolic density profile is assumed. The 
gravitational potentials are initially set to their flat space val- 
ues. After a few iterations, rotation is switched on, and the so- 
lution converges to the stationary and axisymmetric configu- 
ration fixed by the model parameters H c and O. A particular 
approximate solution, obtained from a previous axisymmetric 
one, remains axisymmetric. The sequence of equilibrium mod- 
els is therefore restricted to the subspace of axisymmetric and 
stationary ones which is part of the full configuration space of 
three-dimensional quasi-equilibrium configurations. 

At a certain iteration step J , after convergence is consid- 
ered to be sufficient, a small perturbation 

5v oc e H c r 2 sin 2 6» (cos 2 i/> - sin 2 ^) (42) 

is added to v = In AT, which excites the I — 2, m = ±2 mode. 
Here H c denotes the central value of the log-enthalpy and e is 
a small constant of order 10~ 8 . The three-dimensional gravi- 
tational potentials N and N l respond to this perturbation via 
the field equations and the matter distribution via the first inte- 
gral of motion. The non-axisymmetry of a particular configu- 
ration is conveniently measured by a parameter (3, introduced 
by means of the Fourier expansion of v in ip, 

v (r e , tt/2, ip) = i> (r e , tt/2, 0) + (3 cos 2 V , (43) 

where r c denotes the mean stellar radius in the equatorial 
plane. As mentioned above, C applied to axisymmetric con- 
figurations is contractive in some neighbourhood of the pre- 
viously constructed axisymmetric solution. Axisymmetric per- 
turbations will thus decay exponentially. This may be different 
for the non-axisymmetric perturbation (42), depending on the 
influence of the three-dimensional terms in the field and matter 
equations. The stability of the axisymmetric model is decided 
by inspection of the behaviour of the non-axisymmetry param- 
eter (3 during the continued iteration. Having in mind that we 
operate in the linear regime, we may introduce k, the ampli- 
fication factor of the non-axisymmetry parameter (3 between 
two successive iterations, and so the following three cases can 
be distinguished: 

1 . k < 1 : (3 decreases exponentially and the perturbed config- 
uration converges to the non-perturbed axisymmetric con- 
figuration — the configuration is secularly stable, 

2. k — 1 : (3 does not change during the subsequent iterations 
— the configuration is secularly meta-stable, 

3. k > 1 : (3 grows exponentially and the perturbed configura- 
tion evolves subsequently away from the unstable axisym- 
metric configuration towards a new stable triaxial quasi- 
equilibrium configuration — the configuration is secularly 
unstable. 

To infer the actual stability of a certain configuration, one has 
of course to keep in mind the approximate character of our per- 
turbed equations. In particular, for fully relativistic configura- 
tions, relativistic terms beyond the current approximation level 



of order 1/2-PN which are not included in the present scheme 
will possibly alter the stability against the triaxial secular insta- 
bility in some a priori unknown sense. 

3. Numerical code 

3.1. Pseudo-spectral method 

The numerical implementation of the analytic model described 
in Sect. 2 relies on a pseudo-spectral code which has been de- 
rived from the three-dimensional code presented in Bonazzola 
et al. (1996), but enhanced by the fully three-dimensional treat- 
ment of the shift vector N l . The previous code was itself an 
extension of an originally axisymmetric code which had been 
conceived to construct high precision models of rapidly rotat- 
ing stars in general relativity (Bonazzola et al. 1993) and in 
the following applied to study neutron star models based on 
realistic equations of state (Salgado et al. 1994) and to model 
neutron stars provided with a strong magnetic field (Bocquet et 
al. 1995). We refer to Bonazzola et al. (1993) for a detailed de- 
scription of the numerical techniques which apply essentially 
to the present three-dimensional code as well: 

The various quantities are expanded in Fourier series in tp 
and 8 and in Chebyshev series in r. The elliptic field equations 
are solved after expansion of the actual sources in terms of the 
angular eigenfunction bases of spherical harmonics Y™ (8, ip) 
for A3 and (cos 18, sin 18) for A2. The resulting systems of or- 
dinary differential equations in the radial variable are conve- 
niently solved in coefficient space. The elliptic equations are 
solved exactly in so far as the computational domain covers the 
whole space. This allows to satisfy the exact boundary condi- 
tion of asymptotic flatness at spatial infinity and the proper cal- 
culation of the source terms which fill the entire space. In this 
way, the limitations of approximate boundary conditions such 
as a Robin boundary condition at some finite radius R (York & 
Piran 1982) can be completely overcome. 

Two grids are used to cover the numerical domain. The in- 
ner one embodies the star whereas the surrounding space is 
compactified thanks to a change of the radial variable u = r -1 , 
mapping it onto the finite exterior grid. 

We further want to point out that the present study impres- 
sively demonstrates the ability of a spectral method to exploit 
the particular nature of the problem where only a small num- 
ber of angular modes in the tjj variable is present. As a conse- 
quence, such low a number as N$ = 4 in the interval [ 0, 27T ] is 
sufficient to exactly represent the various quantities up to linear 
order in (3. 

3.2. Vector Poisson equation 

As outlined in Sect. 2.3, a major issue for the computation of 
triaxial quasi-equilibrium configurations at the present level of 
approximation is the calculation of the shift vector N l by so- 
lution of a generalized vector Poisson type equation of the fol- 
lowing form 

AV l + aV l (WiV l ) = S l , (44) 
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where a is a constant number. For — (44) yields an ellip- 
tic equation for V 1 , conducting to a well posed boundary value 
problem. The particular case with a = 1/3, hence 



AV i + | V^ViV') = S l , 



(45) 



is generic part of the York procedure for the solution of the ini- 
tial value problem of general relativity (York 1983). Therefore, 
numerous attempts of a numerical solution of (45) have been 
made in the past. Bowen (1979) has suggested a simplified set 
of equations by solving for the auxiliary variables W l and U 
defined through 



V 1 = W l 



(46) 



With this definition, (45) reduces to two Poisson equations for 
W l and U, namely 



AW 1 =S l , AU = WiW 1 . 



(47) 



This approach has been widely used in analytical and numerical 
work. Note, however, that Evans (1984) was able to construct a 
suitable Green's function in Cartesian coordinates. It incorpo- 
rates the boundary condition V % = at spatial infinity and reads 
explicitly 



7 



(Xi-x'JiXj-x'A 



7 \x — x 



l\2 



Gij (x, x ) 

32n\x — x' 

where Sij is the flat space metric tensor. Gij satisfies 

[ AGij + i V, (V'Gy)] (*-*') - M (3) • 



,(48) 



(49) 



The elliptic equation (45) is hence transformed into the integral 
equation 



V l {x) = ( G i i(x — x') S l (x') d 3 x' . 

JR 3 



In general, calculations have been carried out in Cartesian type 
coordinates which considerably simplify the solution of the el- 
liptic equations (47), because the vector Poisson equation re- 
duces to independent scalar Poisson equations for the individ- 
ual components V x , V v , and V z . Their use allows further to 
circumvent the difficulties associated with the coordinate sin- 
gularities of spherical or cylindrical coordinates when using fi- 
nite difference methods. Nevertheless, in astrophysical appli- 
cations the use of spherical or cylindrical coordinates are of- 
ten much more adapted to the geometry of the problem where, 
for instance, one encounters a matter distribution with compact 
support. Note, however, that because the choice of the local 
vector basis is arbitrary, one might employ Cartesian compo- 
nents that are functions of the spherical coordinates (r, 6, (f>) 
as well. This would allow to adopt Bowen's scheme (47) unal- 
tered. Boundary conditions may, however, take a simpler form 
for the spherical components V r , V e , and V^. 

We propose a new solution method that has been designed 
to fit specifically to problems in a spherical coordinate system 



when the use of spherical vector components is desirable. Be- 
cause our spectral method takes into account the intrinsic reg- 
ularity properties of analytic functions, the first step consists in 
decoupling the (irrotational) scalar from the (divergence-free) 
pure vector part of the involved vector quantities in order to 
obtain a system of ordinary Poisson equations analogous to 
(47). We introduce vector fields V 1 and S l , representing the 
divergence-free part of V 1 and S l respectively, as well as two 
scalar potentials and $. Appropriate boundary conditions 
supplied, one has a unique decomposition 



V 1 = V l + V** , S l = S l + V l $ . 



(51) 



In a first step, we obtain a Poisson equation for <j> by computing 
the divergence of (51) 



A$ = \7iS l . 



(52) 



Equation (44), expressed in terms of the new variables, then 
reads 



AV l + V'((l+a) A* - *) = S l 



(53) 



Taking the divergence of (44), H=((l+a) A* - $) turns out 
to be a harmonic function, thus satisfying AH = 0. We seek 
regular and bounded solutions to the initial equation in R 3 and 
therefore choose H = 0. Consequently, is determined by the 
Poisson equation 



(1- 



A* = $ 



(54) 



Combining (54) with (53), we obtain a vector Poisson equation 
for V 1 which is just (44) applied to divergence-free vector fields 
( ' and.S'. 



(50) AV l = S* 



(55) 



Equations (54) and (55) resemble (47), but the additional con- 
straint ViV 1 =0 considerably simplifies the solution of (55) in 
spherical coordinates when the standard orthonormal frame is 
used. No distinction between contra- and covariant vector com- 
ponents is made in this case. We further drop the tildes on top 
of the vector indices, which had been introduced in Sect. 2.3 to 
distinguish the pseudo-physical components. The explicit form 
of AV l =S i then reads 



A 3 V r 



2V 



2_ fd_ , 1 

r 2 \ 86 tan 8 



V 



(56) 



d<t> 



= S r 



A,V° 



V° 



r 2 sin 2 9 ' 



2 dV r 
' 2 ~df 



(57) 



dVt> 



r 2 sin 9 tan 9 d(j) 
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Fig. 1. iV^ for a rapidly rotating Kerr black hole at a/M = 0.99 where 
m denotes the radius of the horizon. The different curves correspond 
to various values of 9 ranging from 0° to 90°. The asterisks indicate 
the sub-domain boundaries. 



A,V* 



+ 



2 dV r 
1 - 

r 2 sin 2 9 r 2 sin 9 dip 

2 dV e 



(58) 



r 2 sin 9 tan 9 



The solution of (56) to (58) is obtained in the following way: 
combining the divergence-free condition on V % , 

dV r 2V r 1 dV e V 1 &V+ n /IW 

^ + — + -^r + ^— a + ^- a ^T = ' (59) 

or r r oO rta.ii 9 rsmO dip 

with (56) eliminates all components but V r and leads to an or- 
dinary scalar Poisson equation for the auxiliary variable V r = 

rV r , 



A 3 V r = rS r 



(60) 



In order to simplify (57) and (58) we express V e and V* ac- 
cording to 



v e = i0U 
r 08 



dW 



rsin( 



, V 4, 



1 dW 
r~dW 



1 dU 



rsa\9 



,(61) 



where U and W are two auxiliary scalar potentials. Making 
once more use of (59), the original equations (57) and (58) are 
replaced by 



18_ 

r 89 



d 2 U dV r 



dr 2 dr 







r sin 9 dip 



d 



r sin 9 dip 



d 2 U dV r 



dr 2 dr 



+ 



ld_ 

rd9 



A 3 W- 
A 3 W~ 



2 dW 
r dr 

2dW 
r dr 



- S", (62) 



S+. (63) 



An equation that only involves (d r U -V r ) is obtained after 
differentiation of (62), multiplied by r sin 9, with respect to 9 



Fig. 2. Absolute error in corresponding to Fig. 1. 



and of (63), multiplied by r, with respect to <p>. After integration 
over r where the integration constant is set to zero in order to 
make the source term vanish at spatial infinity, one is left with 
an equation that only involves the angular part of the scalar 
Laplacian, namely 



1 



d*_ 

d9 2 



1 d 



1 



d 2 



dip 2 



dU 
dr 



_yr 



= -S r . (64) 



Since V r has already been determined by means of (60), U can 
be computed immediately after solution of (64). In order to fix 
the lacking potential W, an ordinary integration of (63) over ip 
is carried out, and the integration constant is set to zero to com- 
ply with the required vanishing behaviour of the source terms 
of the resulting equation. The final equation for W r , defined by 
W = rW, then becomes 



A 3 W 



f 

Jo 



S 4, 



1 d fd 2 U dV r 



r sin 9' dip \ dr 2 dr 



d9' . (65) 



The system of equations constituted by (60), (64), and (65) is 
equivalent to (55) and significantly simplifies the numerical so- 
lution. Note the absence of S e in the source terms of the final 
equations. This reflects the dependency of the components of 
S i due to the constraint VjS 1 ' = 0. 

Based on the above approach, we have realized a numer- 
ical scheme that solves (44) in a multi domain configuration. 
The computational domain may be extended to spatial infinity 
thanks to a suitable transformation of the radial variable in the 
exterior zone (see Bonazzola et al. (1993) for a previous appli- 
cation of this method). 

Numerical tests based on simple analytic functions revealed 
the well known evanescent error characteristic for spectral 
methods where the numerical errors quickly reach the round- 
off limit of ~ 10~ 14 of the employed machine. It has further 
been applied to solve the general relativistic field and matter 
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equations for numerical models of self-gravitating matter sur- 
rounding a rotating black hole which will be presented else- 
where (Bonazzola & Gourgoulhon, in preparation). For vanish- 
ing matter density, the Kerr vacuum solution is recovered and 
allows a direct comparison with the analytic solution expressed 
in quasi-isotropic coordinates. The coordinate base shift vec- 
tor component near the horizon of a rapidly rotating black 
hole, characterized by a/M = 0.99, is shown in Fig. 1 whereas 
the involved numerical error is illustrated in Fig. 2. The out- 
side of the black hole is covered by three grids where the ex- 
terior compactified one extends to spatial infinity. The actual 
grid resolution is N r = 33 in each zone and Ng = 7, covering 
the interval < 9 < n/2. The numerical error committed on 
nowhere exceeds 10~ 10 with respect to the analytic solu- 
tion and decreases even further to ~ 10~ 13 for higher values 
of N r . Note the very high precision of the numerical model 
in spite of the low number of collocation points, which shows 
the power of spectral methods for problems where the involved 
quantities can be represented by analytic functions. 

3.3. Validation of the code 

The numerical results in the axisymmetric case have been veri- 
fied repeatedly by different tests, based on comparison with an- 
alytic solutions or numerical models provided by other authors 
(Eriguchi et al., in preparation), as well as by exploiting in- 
trinsic error indicators like the general relativistic virial identi- 
ties (Bonazzola 1973; Bonazzola & Gourgoulhon 1994; Gour- 
goulhon & Bonazzola 1994) or analytic properties of involved 
quantities. Typical global errors for rapidly rotating neutron 
star models based on a 7 = 2 polytropic equation of state are 
about 10 -6 . This value has to be confronted with a global error 
of 10~ 14 for non-rotating, spherical configurations. This differ- 
ence stems from a certain deficiency of the spectral approxima- 
tion of non-analytic functions. In the spherical static case, the 
boundary of the inner spherical grid can be chosen to coincide 
with the star surface. As a consequence, all quantities are an- 
alytic in the adjacent sub-domains, and the numerical error is 
dominated by the roundoff error of the employed machine for 
a moderate number of grid points thanks to the exponential de- 
crease of the residual numerical errors as a function of the num- 
ber of grid points. For non-analytic functions, the convergence 
slows down considerably and is worst for discontinuous func- 
tions — this behaviour corresponds to the well known Gibbs 
phenomenon of Fourier series. In the rotating case, the star sur- 
face is located inside the central grid zone. Therefore, matter 
variables like the energy density, whose derivatives are discon- 
tinuous across the star boundary and even become singular for 
7 > 2, are no more analytic functions in this domain and suffer 
hence from the above deterioration. For the current study which 
is only concerned with a polytropic equation of state p = up 1 , 
this signifies that 7 has to be limited to values for which the 
spectral approximation still works reasonably well. In practice, 
7 is bounded by a maximal value 7 max = 3. In particular, in- 
compressible fluid configurations such as the Maclaurin or Ja- 
cobi spheroids are beyond the scope of the present code. Since 



the range of astrophysical interest is still well within in the cur- 
rent limits, these restrictions are not prohibitive. 

Because of the lack of concurrent studies in the relativis- 
tic regime, the validation of the non-axisymmetric part of the 
code is limited to the Newtonian limit where some results 
for compressible and incompressible fluids are available, ei- 
ther obtained by analytic or numerical studies. The different 
tests (Bonazzola et al. 1996) include in particular a confir- 
mation of James' value of the critical adiabatic index, yield- 
ing 7 cri t = 2.238 ± 0.002, which is in perfect agreement 
with the predicted value. A classical result from the theory of 
Newtonian incompressible fluid bodies is the particular value 
T/I^lcrit = 0.1375 at the bifurcation point of the secular bar 
mode instability. Since we cannot verify this result directly be- 
cause of the above restrictions imposed on the equation of state 
(the incompressible case corresponds to 7 — > 00), we studied 
the dependence of T/| W| for increasing 7 and found the suc- 
cessive values to converge to the asymptotic one. A compar- 
ison with results of Ipser & Managan (1981) and Hachisu & 
Eriguchi (1982) in the compressible polytropic case for a fixed 
7 = 3 showed an agreement of better than 0.5% for the crit- 
ical angular velocity f2 K and of about 2% for the values of 
T/\W\ crit (cf. Table 1 of Bonazzola et al. 1996). 

The solution of the 3D-shift vector equation is computed 
following the method presented in Sect. 3.2 and its imple- 
mentation is completely independent from the former scheme 
where only one axisymmetric equation for the only non- 
vanishing shift vector component was solved. The devia- 
tion of for identical models is of order 10 -6 and corre- 
sponds to the global numerical error of a particular configura- 
tion labeled by the model parameters H c and £1 

4. Results for Polytropes 

As mentioned in Sect. 1, earlier investigations in the Newto- 
nian or at most post-Newtonian regime in general made use 
of polytropic equations of state which represent simplified but 
within certain limits reasonable and consistent templates of re- 
alistic equations of state. The critical adiabatic index 7 cr j t is 
defined as the value of 7 for which the critical angular veloc- 
ity O cr ;t coincides with the Keplerian angular velocity Ok, the 
maximum angular velocity supported by a rapidly rotating star 
without mass shedding from the equator. For 7 < 7 cr it, the star 
thus cannot rotate rapidly enough for the secular instability to 
develop. If 7 > 7 cr it> there always exists a certain f2 cr it < Ok- 
Provided some mechanism that fuels the star with rotational 
kinetic energy, it may be affected by the secular instability. 

Newtonian polytropic stars obey a scaling law and as a con- 
sequence, 7 cr ; t is a global constant that does not depend on 
other model parameters. For relativistic configurations, effects 
of general relativity are likely to influence the bar mode in- 
stability. Therefore, 7 cr ;t will be a function of the strength of 
the relativistic fields which we measure by the value N c of the 
lapse function at the centre of the star. Fig. 3 shows this de- 
pendence for the case of the 2D-shift vector equation of our 
previous investigation as well as the modified graph for the 
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present approach which includes the solution of the 3D-shift 
vector equation. 7 cr ;t increases only slightly in the first case. 
One even notes a tiny decrease in the weak field regime. For 
the present approach, the situation changes quite dramatically. 
In the Newtonian regime, both graphs converge to James' clas- 
sical value, which, of course, is required by consistency: for 
decreasing field strength the shift vector contribution dies out 
more rapidly than that of the lapse function, so that in both 
cases the asymptotic space-time line element is given by (11). 
In the relativistic domain, however, the shift vector shows to 
inhibit the secular instability more and more efficiently. The 
calculations are not continued beyond 7 = 3 where the steepen- 
ing gradients of the matter variables severely affect the spectral 
approximation. 

The question which arises immediately addresses the origin 
of the very different behaviour for the both approximations. Let 
us recall the classical limit of the first integral of motion 

h + U - \ O V = const . (66) 

The evolution of a perturbed configuration is driven by the 
mutual influence of the gravitational potential U and the en- 
thalpy h which are the only three-dimensional quantities in this 
case. The centrifugal potential is axisymmetric and remains un- 
changed. Since relativistic effects tend to increase the overall 
gravitational forces, one might expect that relativistic correc- 
tions to U inhibit the symmetry breaking. However, the lower 
graph in Fig. 3, which illustrates the results of our previous in- 
vestigation, shows an only slight growth of 7 cr ;t and, moreover, 
even a destabilizing effect in the weakly relativistic regime. In 
the general relativistic case, (66) is replaced by 

H + v - In r = const . (67) 

The crucial point is that now the "centrifugal" potential Uq = 
— In r depends on the azimuthal angle as well. This means, the 
angular modulation of Uq — v is superimposed by that of Uq, 
the contribution from the latter being a purely relativistic ef- 
fect without Newtonian counterpart. The variation of Uq with 
respect to TV may then enhance or diminish the triaxial pertur- 
bation of the effective relativistic potential, depending on its 
sign and, via the first integral of motion, prevent or favour the 
triaxial deformation of the fluid body. The relativistic centrifu- 
gal potential Uq according to (37) reads 

£/fi = -lnr = ±ln(l-*7|), (68) 

where U^, the physical fluid velocity in direction as measured 
by the local Eulerian observer O , specified in (35), has the 
form 

In the 2D-shift vector case, the only three-dimensional quantity 
in (69) is the lapse function N. The variation SUq of Uq with 
respect to N adds to 6Uq with a positive sign: the triaxial per- 
turbation of the relativistic potential is enhanced. This explains 




1 0.8 0.6 0.4 

Fig. 3. Critical adiabatic index 7 cr it as function of N c , the lapse func- 
tion at the centre of the star. The upper curve is obtained with solution 
of the full 3D-shift vector equation whereas the lower one illustrates 
the dependence for the axisymmetric computation of the shift vector. 
The intermediate curve shows the analogous curve for the first case 
with the axisymmetrized centrifugal potential. 

the destabilizing effect in the weak field regime before the over- 
all growth of the relativistic forces related to Uq dominates in 
the strong field region. The computation of the shift vector in 
a three-dimensional fashion leads to an additional variation of 
the term (ft — N^) which contributes to 5Uq with a negative 
sign. It competes with the contribution of the lapse function N 
to SUq, acting in the opposite sense. It turns out that the com- 
bined effect of both contributions is dominated by the stabiliz- 
ing effect of the shift vector, which corresponds to the upper 
curve in Fig. 3. From the Newtonian limit on, the symmetry 
breaking is increasingly suppressed. The exclusive influence of 
the relativistic gravitational potential Uq can be studied by av- 
eraging Uq over ip. This corresponds to the Newtonian situa- 
tion, and the resulting curve for the critical adiabatic index in 
this case has been added to Fig. 3. It is located between the two 
other curves, obtained with a three-dimensional Uq, and vali- 
dates our above reasoning. We have confirmed numerically that 
the fluid velocity U^, appearing in the relativistic centrifugal 
potential Uq, is the only quantity where the three-dimensional 
character of the shift vector actually affects the bar mode insta- 
bility. We have further checked the robustness of our results by 
adding further three-dimensional terms to the field equations 
due to an additional extra-diagonal element in the metric ten- 
sor or to treating A and B as "pseudo" 3D-variables. No no- 
ticeable modification of the results has been observed. These 
numerical tests give us confidence in the analytic approxima- 
tion and indicate that an additional 3D-treatment of the tensor 
part of the space-time metric tensor is unlikely to modify the 
present results significantly. 

As concerns Wilson's approach, we have carried out the 
above calculations by imposing A = B which mimics the 
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"conformally flat condition" except that we do not solve for 
the fully three-dimensional conformal factor but rather adopt 
A 2 N~ 2 as conformal factor where A remains axisymmetric. 
The global numerical error of the initial axisymmetric config- 
urations as measured by means of the general relativistic virial 
theorems is of order 1CU 4 which has to be compared with val- 
ues of 10 -6 for solution of the exact equations. The results 
for the critical adiabatic index 7 cr ; t appear to be rather insen- 
sitive to this simplification. The actual values coincide with 
those by solution of the original equations within a relative er- 
ror of about 1CP 3 and confirm the validity of Wilson's sim- 
plified scheme for quasi-equilibrium configurations at least for 
this particular application. 

5. Conclusion 

The present work has studied the viscosity driven secular in- 
stability of rapidly rotating stars in general relativity by solu- 
tion of an approximate set of field equations related to Wil- 
son's scheme for quasi-equilibrium configurations. It has been 
revised to admit a non-conformal metric tensor which allowed 
us to recover the exact equations for stationary and axisym- 
metric configurations in contrast with the simplifying "confor- 
mally flat condition". The main improvement with respect to 
our previous study consists in the fully three-dimensional treat- 
ment of the shift vector. The numerical scheme has been ap- 
plied to configurations built on a polytropic equation of state. 
The shift vector shows to strongly enhance the stability of rela- 
tivistic configurations compared to the contribution of the lapse 
function only. The weak variation of the critical adiabatic index 
7 cr it in the previous investigation can be explained by a par- 
tial cancellation of certain non-axisymmetric terms in the first 
integral of motion. Comparative calculations adopting the sim- 
plifying "conformally flat condition" essentially yield the same 
results. This is intelligible, having in mind that even for maxi- 
mum rotation neutron star models the deviation of the geometry 
of curved space from conformal flatness is of order 10~ 3 . 

Nevertheless, the symmetry breaking is still possible for a 
mass of M ns = 1.4M and an adiabatic index of 7 = 2.5. A fu- 
ture investigation has to address the question, if realistic models 
of neutron star matter still admit the symmetry breaking for as- 
trophysically relevant masses, as it was the case in our previous 
investigation with an axisymmetric shift vector. 
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